/*
 * MineralKinetics.cpp
 *
 *  Created on: 21 Jul 2011
 *      Author: Allan
 */

#include "MineralKinetics.h"

// C++ includes
#include <cmath>
using namespace std;

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include <Core/Multiphase.h>
#include <Numerics/AuxiliaryFunctions.h>
#include <Utils/Assert.h>

MineralCatalyst::MineralCatalyst()
{}
	
MineralCatalyst::MineralCatalyst(const string& name, const string& quantity, double power)
: name(name), quantity(quantity), power(power) 
{}

bool 
MineralCatalyst::operator==(const MineralCatalyst& mineral_catalyst) const
{
	if(name != mineral_catalyst.name)
		return false;
	
	if(index != mineral_catalyst.index)
		return false;
	
	if(quantity != mineral_catalyst.quantity)
		return false;
	
	if(!AlmostEqual(power, mineral_catalyst.power))
		return false;
	
	return true; 
}
	
bool
MineralCatalyst::operator!=(const MineralCatalyst& mineral_catalyst) const
{
	return !(*this == mineral_catalyst);
}

MineralMechanism::MineralMechanism()
{}

MineralMechanism::MineralMechanism(double k, double Ea)
: k(k), Ea(Ea)
{}

MineralMechanism::MineralMechanism(double k, double Ea, const MineralCatalyst& catalyst)
: k(k), Ea(Ea), catalysts(1, catalyst)
{}

MineralMechanism::MineralMechanism(double k, double Ea, const vector<MineralCatalyst>& catalysts)
: k(k), Ea(Ea), catalysts(catalysts)
{}

bool
MineralMechanism::operator==(const MineralMechanism& mineral_mechanism) const
{
	if(!AlmostEqual(k, mineral_mechanism.k))
		return false;
	
	if(!AlmostEqual(Ea, mineral_mechanism.Ea))
		return false;
	
	if(catalysts != mineral_mechanism.catalysts)
		return false;
	
	return true;
}

bool
MineralMechanism::operator!=(const MineralMechanism& mineral_mechanism) const
{
	return !(*this == mineral_mechanism);
}

MineralKineticsInfo::MineralKineticsInfo() 
{}

MineralKineticsInfo::MineralKineticsInfo(const string& mineral) :
mineral(mineral)
{}

void
MineralKineticsInfo::SetMineralName(const string& mineral)
{
	this->mineral = mineral;
}

void
MineralKineticsInfo::SetSpecificSurfaceArea(double specific_surface_area)
{
	this->specific_surface_area = specific_surface_area;
}

void
MineralKineticsInfo::AddMineralMechanism(const MineralMechanism& mechanism)
{
	mechanisms.push_back(mechanism);
}

void
MineralKineticsInfo::SetMineralMechanisms(const vector<MineralMechanism>& mechanisms)
{
	this->mechanisms = mechanisms;
}

const string
MineralKineticsInfo::GetMineralName() const 
{
	return mineral;
}

const double
MineralKineticsInfo::GetSpecificSurfaceArea() const 
{
	return specific_surface_area; 
}

const vector<MineralMechanism>&
MineralKineticsInfo::GetMineralMechanisms() const 
{
	return mechanisms;
}

MineralKinetics::MineralKinetics() :
imineral(-1), specific_surface_area(0.0)
{}

MineralKinetics::MineralKinetics(const MineralKineticsInfo& info, const Multiphase& multiphase) :
mineral(info.GetMineralName()), mechanisms(info.GetMineralMechanisms()), 
imineral(multiphase[mineral]), specific_surface_area(info.GetSpecificSurfaceArea())
{
	// Set the indexes of the catalysts species in the mineral mechanisms
	BOOST_FOREACH(MineralMechanism& mechanism, mechanisms)
		BOOST_FOREACH(MineralCatalyst& catalyst, mechanism.catalysts)
			catalyst.index = multiphase[catalyst.name];
}

MineralKinetics::~MineralKinetics()
{}

const double
MineralKinetics::KineticRate(const Reaction& reaction, double T, double P, const VectorXd& n, const VectorXd& a) const
{
	// The temperature in Kelvin
	const double Tk = T + 273.15;
	
	// The universal gas constant in kJ/(mol.K)
	const double R = 8.314472E-3;
	
	// The equilibrium constant of the mineral reaction
	const double K = reaction.EquilibriumConst(T, P);
	
	// The reaction quotient of the mineral reaction
	const double Q = reaction.ReactionQuotient(a);
	
	// The reactive surface area of the mineral, which is given by the product of the mineral specific surface area and its current number of moles
	const double A = specific_surface_area * n[imineral];
	
	// The effective rate constant taking into account all mechanisms
	double keff = 0.0;
	
	// Iterate over all the mechanisms
	BOOST_FOREACH(const MineralMechanism& mechanism, mechanisms)
	{
		// The kinetic rate constant for the current mechanism at 25 C
		const double k25 = mechanism.k;
		
		// The Arrhenius activation energy for the current mechanism
		const double Ea = mechanism.Ea;
		
		// Calculate the rate constant for the current mechanism
		double k = k25 * exp(-Ea/R * (1.0/Tk - 1.0/298.15));
		
		// Update the rate constant for the current mechanism with the catalyst contributions
		BOOST_FOREACH(const MineralCatalyst& catalyst, mechanism.catalysts)
			k *= pow(a[catalyst.index], catalyst.power);
		
		// Update the effective rate constant taking into account all mechanisms
		keff += k;
	}
	
	// Calculate the rate of dissolution/precipitation of the mineral (positive: dissolution, negative: precipitation)
	return keff * A * (1.0 - Q/K);
}

const VectorXd
MineralKinetics::DerivativeKineticRate(const Reaction& reaction, double T, double P, const VectorXd& n, const VectorXd& a) const
{
	return VectorXd::Zero(n.rows());
}
